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1 Introduction 

We consider high-order central approximations for solutions of one-dimensional 
Hamilton- Jacobi (HJ) equations of the form 

+ = °> xeR ’ (1) 

subject to the initial data <j>{x,t = 0) = &>(x). Solutions for (1) with smooth 
initial data typically remain continuous but develop discontinuous derivatives 
in finite time. Such solutions are not unique; the physically relevant solution 
is known as the viscosity solution (see [1, 3, 4, 5, 8, 15] and the references 

therein). _ . 

Various numerical methods were proposed in order to approximate the so- 
lutions of (1). Examples for such methods are the high-order Godunov- type 
schemes that were introduced in [20, 21], and were based on an Essentially 
Non-Oscillatory (ENO) reconstruction step [7] that was evolved in time with 
a first-order monotone flux. The least dissipative monotone flux, the Go- 
dunov flux, requires solving Riemann problems at cell interfaces. A fifth-order 
Weighted ENO (WENO) scheme, based on [10, 18], was introduced by Jiang 
and Peng [9]. 

Recently, Lin and Tadmor introduced in [16, 17] central schemes for ap 
proximating solutions of the HJ equation. These schemes are based on the 
Nessyahu- Tadmor scheme for approximating solutions of hyperbolic conser- 
vation laws [19]. Unlike upwind schemes, central schemes do not require Rie- 
mann solvers, which makes them attractive for solving systems of equations 
and for multi-dimensional problems. A second-order semi-discrete version 
of these schemes was introduced by Kurganov and Tadmor in [12]. While 
less dissipative, the semi-discrete scheme requires the estimation of the lo- 
cal speed of propagation, which is computationally intensive in particular m 
multi-dimensional problems. In a later work [11], the numerical viscosity was 
further reduced by computing more precise information about local speed of 
propagation. To address the problem of schemes that are too computation- 
ally intensive, we introduced in [2] efficient first- and second-order central 
schemes for approximating the solutions of multi-dimensional versions of (1). 
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Unlike the previous attempts, our schemes in [2] scale well with increasing 

dimension. / 

In this paper we derive fully-discrete Central WENO (CWENO) schemes 
for approximating solutions of (1), which combine our previous works [2, 
13 14]. We introduce third- and fifth-order accurate schemes, which are the 
first central schemes for the HJ equations of order higher than two. The 
core ingredient in the derivation of our schemes is a high-order CWENO 
reconstructions in space. 

Acknowledgment: We would like to thank Volker Elling for helpful discus- 
sions. 


2 CWENO Schemes for HJ Equations 


We are interested in approximating solutions of (1) subject to the initial data 
4>(x t = 0) = <j>o(x). For simplicity we assume a uniform grid grid in space 
and time with mesh spacings, h := Ac and At. We denote the grid points by 
x = iAx t n = nAt, and the fixed mesh ratio by A = At/ Ax. Let denote 
the approximate value of 4>{x u t n ), and (p*)? denote the approximate value of 
the derivative 4> x (Xj,t n )- We define A + ip^ := y^+1 — Vi > ^ Vi Vi Vi - 1 


and A°(fii '■= V?+i ~ V?-i- . . . . T 

We assume that the approximate solution at time t", ^ given. In 

order to approximate the solution at the next time step t , Vi > start 
by reconstructing a continuous piecewise-polynomial from the data, <f/. and 
sample it at the half-integer points, {x i+ i/ 2 }, in order to obtain the point- 
values of the interpolant at these points </?? + i/ 2 as wel1 35 the denvatlve > 


<p' i+1/2 . We then evolve from time 


t n to time t n+l according to (1), 


V 


£ Tt “J— 1 

(x i+ r,C +1 ) = (x i+ l,f n ) - [ n H (vr (*<+*•*)) 


dt. 


( 2 ) 


This evolution is done at the half-integer grid points where the reconstruc- 
tion is smooth (as long as the CFL condition A \H' (<^)| < 1/2 is satisfied). 
Finally, in order to return to the original grid, we project V^+i/2 ' Dack onto 

the integer grid points {x,} to end up with 

Since the evolution step (2) is done at points where the solution is smooth, 
we can approximate the time integral at the RHS of (2) using a sufficiently 
accurate quadrature rule. For example, for a third- and fourth-order method, 
this integral can be replaced by a Simpson s quadrature, 

’ H (vx (x 1+ r,t)) * * f [* (V- (*<+*■*")) 

+4 H (x j+ l,t" + 5)) +H (<Pr (x i+ I,t" +1 ))] - 


(3) 
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The intermediate values of the derivative in time, (x l+1 / 2 , t n+1 ‘ ), and 
V (x + 1/2 t n+1 ), which are required in the quadrature (3), can be predicted 
using a Taylor expansion or with a Runge-Kutta (RK) method. For details 
we refer the reader to [13, 19] and the references therein. 

The remaining ingredient is the piecewise- polynomial reconstruction in 
space A careful study of the above procedure reveals that there are actually 
three different quantities that should be recovered in every time step. First, 
given <pi at time t n we need to reconstruct the point- values at the half-integer 
grid points, <^ +1/2 , at the same time t". This is the first term on the RHS of 
(2) The second term on the RHS of (2) requires evaluating the Hamiltonian 
H at the derivative y>' +1/2 . Hence, the second quantity we should recover 
is 2 from tpo Finally, the predictor step that provides the values at 

the quadrature nodes in (3), require us to estimate y?' +1/2 from V*+ 1/2 at 
every step of the RK method. In the next two sections we will focus on the 
reconstruction of these three quantities, first for a third-order method an 
then for a fifth-order method. n+1 . 

The projection from <^7+1/2 onto t!ie ori 8 inal S rid P omts to ls 

accomplished using the same reconstruction used to approximate v”+i/2 ^ ronl 


2.1 A Third-Order Scheme 

Following the above procedure, a third-order scheme can be generated by 
combining a third-order accurate ODE solver in time with a sufficiently 
high-order reconstruction in space. Here we present fourth-order CWENO 
reconstructions of the point values of ¥h + 1/2 ar >d its derivative Vh+i/2' 


The reconstruction of ifii+ 1/2 from <fi- 

In order to obtain a fourth-order reconstruction of ip i+ 1/2 we will write a con- 
vex combination of two quadratic polynomials, tp® constructed on ^ a stencil 
which is left- biased with respect to x i+1 / 2 , and the right-biased 1 p + , 

<^L 21 (x) =Vi + l (A + V>i) (* - xi) + ^ (* - **) (* - x ’+i) + 0 ( h3 ) ’ 

$ (x) = <Pi + \ (* + <Pi) (x - Xi ) + ^2 (A + A + <Pi) ( x - x ’) < x - + ° ^ ■ 

An evaluation of these approximations at {x i+ i } reads 


( x »+i) = ¥’+ ( x ’+s) _ g (3v?i+6<*>i+l <fi+2)- 

A straightforward computation shows that 

(**j) + ^+(*i+j) = + ° ( hA ) ■ 
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The fourth-order WENO estimate of <*p l+l /2 is therefore given by the convex 
combination 



where the weights satisfy ttn +1/ , 2 + w ?+i/2 ^ w t+i/2 — smoot ^ 

regions we would like to satisfy w i ~ ~ \ to attain an O [h ) eiror, 

while when the stencil {xi_i, x t , x i+1 , x i+2 } supporting <p w (x, + i ) contains 
a discontinuity, the weight of the more oscillatory polynomial should vanish. 
Following [10, 18], we meet these requirements by setting 


<k = 


*+5 




l 

i+h 



(4) 


where k,l £ {+,-} (k and l will range over a larger space of symbols when 
we use more interpolants). The constants c* = 1/2 and are independent of 
the grid-point. We choose e as 1(T 6 to prevents the denominator in (4) from 
vanishing, and set p = 2 (see [10]). The smoothness measures Sf should be 
large when ip is nearly singular. Following the standard practice with W ENO- 
type schemes [10], we take Sf 1 to be the sum of the L -norms of the first and 
second derivatives on the stencil supporting . If we approximate the first 
derivative at Xj+ 1/2 by j- t the second derivative by -^A + A <A+ 1 / 2 > 
and define the smoothness measure 


S i+ 1 [r,s\ 


= *E 






+ h 


E 

Jf=r+1 


{b 


A + A~ 


'fi+i+i 


)’■ 


(5) 


then for the fourth-order interpolation of <p w we have S i+1/2 - 

$+ 1/2 [-1,0] and S ?+\/2 = 5 *'+V 2 [M* 


The reconstruction of S^+ 1/2 f rom Vi- 

To obtain a fourth-order estimate of the derivative 1 / 2 ) from we 

start from the cubic interpolants 

y?!? (x) = ^ ^ ( x — x % ) + ^2 ^0 ( x “ “ Xi +i) 

H — (A~A + A“<Pi) (x - Xi) (x - x i+ 1 ) (x - x i _ 1 ) + O ( h 4 ) , 
6h 6 

p { l ] (x) = Pi + i (4 + *») (x-Xi)-f (4 + 4 + *>i) (x - X,) (x - X i+1 ) 

+-J— (A + A + A + ‘f i ) (x-x i )(x-x i+ l)(x-x ^+ 2) + 0(/l‘ , ) . 

6 h 6 

[31 

Differentiating <p± at x i+ i 
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>' |3! - = — {ifi- 2 - 3^_1 - 21^ + 23 ^, +i) , 


-.>+5 


24/i 


Again, 


i = jrrr (—23 (pi 4 21<^i+i 4 3^1+2 ¥^+ 3 ) • 
+ >*+2 24/i 


+ 2^S+i " ^+i + ° ^ ' 


and a fourth-order WENO reconstruction of <p' (^i+i) is 

*£1/2 = w 7 + 1 ^- 3 !+ i + u, «+i ^+!!+ i 

where the weights are of the form (4) with = 1/2 and S i+ i/ 2 “ 
^1+1/2 [”2,0] and S+ +1 / 2 — 1/2 [0 7 2] . 

The reconstruction of ^^1/2 fr° m ¥>141/2- 

Repeating the above procedure, this time with three quadratic interpolants 

V [ - ] (*) = ‘Pi+i + 1 (A~<Pi+ 1) (* ~ x i+0 

+^2 ( x " Xi +0 ( x -^+f) +°( h3 )’ 

vf ( x ) = ¥>«+§ + ^ (a°v? i+ i) (x - *i+i) 

{ A+A ~*i + i) ( x - **-0 ( x - x *+i ) + 0 1 

<£+ ( x ) = Vi+i + £ (^’Vi+i) ( x “ x *+i) 

(^ + ^+0 ( x - x ^+0 ( x - x -+§) + 0 (/j3) ’ 


results with 


where 


V^i , + V 2 ) t A + o (/i 4 ) 

6^-,i+i 3*»M+i 6 +- t+ i * + 2 


rm 1 ~r' ro1 1 , 

0'1‘i+i = 2*^-1 " 4 ^-i + 3v ’*+i ) ’ ^<M+J = 2/4’+§ " ^4’ 

1 




+ 4H 


— (-3^ i+ i +4v? i+ | -?<+$)• 


The fourth-order WENO estimate of 4- +1 / 2 is 

= »r+ j^1 +j + <§CU + 

where the weights re are of the form (4) with c = c + = 1/6, c° = 2/3, and 
the oscillatory indicators S~ +1 ^ 2 = Si+1/2 [ _ 2, -1], S i+1/ / 2 — Si+i/2 
and S+. 1/2 = Si+1/2 [0, 1]. 
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2.2 A Fifth-Order Scheme 

Once again, similarly to the third-order scheme, we need to reconstruct the 
point-values of ip and <p'. We start with the reconstruction of ^.+ 1/2 and 
V^+i/2 from We write sixth -° rder interpolants as a convex combination 

of cubic interpolants, (x) and (x) introduced above and 

p [ o ] (x) = Vi + l (A + <Pi) (x - Xi) + (A + A~p t ) (x - Xi) (x - x 1+1 ) 

+ -L M + 4" A*Vi) (x - Xi) (x - x i+ i) (x - x i+2 ) + O ( h 4 ) . 

6n J 


In this case 


,P] 


16 




, 5 [3] 

+ o^0,i+l 


+ 


n [3] 


16^+d+i 


,+0(h 6 ), 


where 


¥> 13 ' , = 2 - 5v?i-i + 15^ + 5y> i+ i), 

»* + 2 10 

d? 1 . , = A(— i + 9 <fi + 9^ i+ i - Pi+ 2 ), 

u ^+2 lO 

< 2 ^ . , 1 = -L. + I5<pi+i - 5<p i+ 2 + ^+ 3 )- 

1 j 1 1 2 lb 


In a similar way, 


9 '[3] 

-§o <p -.‘+* 



/[3] 

0 ,t+ 5 


9 /[3] 

80 ^+’*+5 


^i+ 1/2 + ^ (^ 6 ) ’ 


where 

, ! = T7r(^i-2 - 3y?i_i - 21^ + 23^ l+ i), 

“>*+3 24 n 

I = ” C ^^ >i + _ 

^l 3 ] = — -- ( — 23^i -f 21^+1 + — ^ 1 + 3 ) 

+i*+2 24/i 

The sixth-order WENO estimates for <fii +\/2 and Pi+ 1/2 are 


i<P 


[3] 


[61 

d+l=^, 2 , 

'(6| ^[3] 

V’i+i - w i+^-,i- 


+ fctf 


,131 

0,i- 

~'[3] 


i+i^0,i+± 


+ <+ 


[3] 


+ < 1^1 +<i^'i 3 i + i’ 


where the weights for <p are given by (4), with c_ = c+ = 3/16, cq = 5/8 and 
the oscillatory indicators are S~ +1 ^ 2 = Sj+i/2 [ 2, 0], S' t n + 1 / 2 = ^b+ 1/2 Ef] 
and 5+_j = S 1+1 / 2 [0,2]. The negative weights for <p' require special treat- 
ment (see [22] for details). Following [22] we split the positive and negative 
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weights in the following way: first, we set 7- = 7± = 9 / 40 > 7° - 49/40 and 
j ~ — y + = 9/80, 7+ = 49/20. Then, For fc, Z € {-,0, +}, set a± = Efc 7± so 
that similarly to (4), 


a 


i.i+5 




P 


and 


tk 

W i+i 


^+ 


+ ,t+i 




E 


i a +,t+i 


Ei a i 




Because <y?[+ 1//2 an<4 ^i+1/2 are defined on the same stencils ’ the ^ use the 
same smoothness measures 4-1/2- 

All that is left is the reconstruction of </ i+ i/ 2 from V’i+i/2- In ,hls case 
a sixth-order approximation to <p' i+1/2 requires a weighted sum of four cubic 
interpolants. This reconstruction is similar to the previous ones. We skip the 
details and summarize the result: 


'[ 6 ] 


*V+ 


W . 


— /[3] 


i+ 


-/[ 3) 

+ 0 -, 


t+i 


+<*$».«+* +w r+* 


, 5 '( 3 1 

0+,1+i 


V 5 , 


,[3] 

+ ,»+*’ 


where 


^l 3 i + i = + 9<(5 *-§ _ 18 ^-| + 11< ^*+§ t 

i ” 6i ^-* + 3l O+i + 

<^o+t + i = + 6 ^+§ - V’.+ f)’ 

= + 18l,5< +§ ~ 9l,p< +f + 2< < c ’i+|)- 

Here, c_ = c+ = 1/20, Co— = Co+ — 9/20, S i+1 ^ 2 = ^t+i/2 — 1], ‘-’,+1/2 — 

S i+ 1/2 [-2,0], S®+ /2 = 5 <+1/a [-1, 1] and 5+ 1/2 = $<+1/2 [0,2], 


3 Numerical Examples 

In all our numerical simulations, the ODE solvers we use are the non-linear 
fourth-order Strong- Stability Preserving Runge-Kutta (SSP-RK) methods of 
[ 6 ], 

We start by testing the accuracy of our new CWENO methods when 
approximating the solution of the linear advection equation, <p t + <p x = 0. 
The initial data is taken as <p(x, 0) = sin 4 +1), the mesh ratio A = 0.9 and 
the time T = 4. The results obtained with the fifth-order method of §2.2 are 
shown in Table 1. 
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Table 1. Error and convergence rate for linear advection with initial condition 
ip(x, 0) = sin 4 (7rx) 


N L\ error 

L\ order 

50 5.03 x 10“ 2 

- 

100 8.36 x 10“ 5 

9.23 

200 2.56 x 10“ 6 

5.03 

400 8.24 x 10' 8 

4.96 

800 2.99 x 10“ 9 

00 


Next, we test the CWENO methods with two nonlinear Hamiltonians, 
a convex Hamiltonian 9 ^ + 5 (‘px + l ) 2 = 0 and a non-convex Hamilto- 
nian Ifit - cos {tp x + 1) = 0 . The interval is [ 0 , 2 ], the boundary conditions 
are periodic and the initial conditions for both Hamiltonians are taken as 
<p(x,0) = — cos (■kx). The exact solution to both problems is smooth until 
t « 1 / 7 r 2 , after which a singularity forms. A second singularity forms m the 

non-convex H example at t w 1.29 /tr 2 . 

The results of the accuracy test with the fifth-order method are shown in 
Table 2 , and the solution at time T = 1-5/tt is plotted in Figure 1. Following 
[ 9 ] the errors in Table 2 after the formation of the singularity are computed 
at a distance of 0.1 away from any singularities. 


Table 2. Lx 

Hamiltonians. 


Error and convergence rate estimates for convex and non-convex 
top: T = 0.5/tt 2 , bottom: T = 1.5/tt 2 . A = 0.3 


N 

convex 
L\ error 

convex 
L\ order 

non-convex 
L\ error 

non-convex 
L\ order 

50 6.35 x 10' 6 
100 1.62 x 10 -7 
200 5.72 x 10' 9 
400 2.73 x 10 _1 ° 
800 1.45 x 10' 11 

5.30 

4.82 

4.39 

4.23 

4.17 x 10' 5 
1.49 x 10' 6 

4.19 x 10' 8 
1.34 x 10' 8 

4.20 x 10 -8 

4,81 

5.15 

4.97 

4.99 

N 

convex 
L\ error 

convex 
L\ order 

non-convex 
L\ error 

non-convex 
L\ order 


a 


0 cc w in -5 



100 1.03 x 10“ 5 
200 9.68 x 10“ 8 
400 6.20 x 10' 10 
snn 1 qo x 10 _U 


4.37 
6.73 
7.29 
?; ns 


7.80 x 10 -7 
1.70 x 10 -8 
5.02 x 10' 10 
1 71 x 10 _u 


5.03 
5.52 
5.08 
i Sfi 
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Fig. 1. left: Convex Hamiltonian right: non-convex Hamiltonian at T — com- 
pared with the exact solution, N = 100. 
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